Draft version August 10, 2009 

Preprint typeset using I^T^X style cmulatcapj v. 02/09/03 



DYNAMICAL EVOLUTION OF THIN DISPERSION-DOMINATED PLANETESIMAL DISKS 

Roman R. Rafikov 1,2 & Zachary S. Slepian 1 

Draft version August 10, 2009 



o 
o 

(N 
OX)! 

< 

o 



PL, 
Ph 1 

6 

o3 



> 

oo 
On 

m 

00 

o 

On 
O 



13 



ABSTRACT 

We study the dynamics of a vertically thin, dispersion-dominated disk of planetesimals with eccentric- 
ities e and inclinations i (normalized in Hill units) satisfying e> 1, i< e~ 2 <C 1. This situation may 
be typical for e.g. a population of protoplanetary cores in the end of the oligarchic phase of planet 
formation. In this regime of orbital parameters planetesimal scattering has an anisotropic character 
and strongly differs from scattering in thick (i ~ e) disks. We derive analytical expressions for the 
planetesimal scattering coefficients and compare them with numerical calculations. We find significant 
discrepancies in the inclination scattering coefficients obtained by the two approaches and ascribe this 
difference to the effects not accounted for in the analytical calculation: multiple scattering events 
(temporary captures, which may be relevant for the production of distant planetary satellites outside 
the Hill sphere) and distant interaction of planetesimals prior to their close encounter. Our calcula- 
tions show that the inclination of a thin, dispersion-dominated planetesimal disk grows exponentially 
on a very short time scale implying that (1) such disks must be very short-lived and (2) planetesimal 
accretion in this dynamical phase is insignificant. Our results are also applicable to the dynamics of 
shear-dominated disks switching to the dispersion-dominated regime. 

Subject headings: 



1. INTRODUCTION . 

Terrestrial planets are thought to be formed by ag- 
glomeration of a large number of primitive rocky or icy 
bodies known as planetesimals (Safronov 1972). While 
the origin of planetesimals themselves is still a rather un- 
certain issue (Youdin 2008) the process of their collisional 
agglomeration has been extensively explored (Wetherill 
& Stewart 1989, 1993; Kenyon & Luu 1998; Kenyon & 
Bromley 2004, 2009). Gravitationally induced bending of 
the trajectories of interacting bodies called gravitational 
focusing (Safronov 1972) is known to play a very im- 
portant role in speeding up the agglomeration process. 
The degree to which the planetesimal collision rate is 
amplified by focusing depends sensitively on the velocity 
dispersion of the planetesimals: the lower is the relative 
velocity between the interacting bodies the higher are 
the gravitational focusing and the collision cross-section. 
Thus, understanding the accretion history of planetesi- 
mals is impossible without understanding their dynami- 
cal evolution. 

Evolution of planetesimal velocities is driven mainly 
by their mutual gravitational interaction. A convenient 
way to characterize the shape of planetesimal orbits and 
their interaction is via the so-called eccentricity and in- 
clination vectors e and i defined as (Ida 1990) 

e = (e x , e y ) = (e cos r, e suit), 

i = (i x , i y ) = (i cosw, i sinw), (1) 

where e, i, r and to are, respectively, the eccentricity, in- 
clination, and horizontal and vertical phases of the plan- 
etesimal. Scattering of two low-mass planetesimals de- 
pends only on their relative eccentricity and inclination 
vectors e r = ei — e2 and i r = L — 12 (the so-called Hill 
approximation, see Henon & Petit 1986). 
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There are two important asymptotic regimes of plan- 
etesimal interaction: shear- dominated and dispersion- 
dominated. The former is realized when the random com- 
ponent of planetesimal velocity, determined by its eccen- 
tricity and inclination, is small compared to the Hill ve- 
locity vh = 0,Rh- Here fi is the local angular frequency 
in the disk, Rh = a(fxi + A^) 1 ^ 3 is the Hill radius, de- 
termined by the distance a to the central object and the 
masses of the interacting bodies ui\ and ni2 relative to 
the central mass M±: /ij = rrii/M*, i — 1,2. Introducing 
scaled relative eccentricity 3 e r and inclination i r vectors 
of the interacting bodies as 4 e r = e r /(/ix + A^) 1 / 3 and 
i r = i r /(/ii + z^) 1 ^ 3 , one can rewrite the condition for 
the shear- dominated regime as 



el+~il< 1. 



(2) 



The relative speed of a pair of interacting bodies in this 
regime is set mainly by the Keplerian shear. 

Dispersion-dominated regime of planetesimal interac- 
tion is realized when 



el + il > 1. 



(3) 



In this case the relative velocity of the planetesimals 
is determined mainly by their random epicyclic mo- 
tion while Keplerian shear plays only a minor role. 
This makes possible analytical treatment of planetesi- 
mal dynamics (Ida 1990; Ida & Makino 1992; Tanaka k 
Ida 1996,1997; Stewart & Ida 2000), which until now 
has been concentrated on the case when i ~ e, so 
that the random velocity distribution of planetesimals 
is roughly isotropic. This assumption is very natural in 

3 In this paper all quantities with a tilde are assumed to be 
scaled by the Hill factor (/xi + ^t2) 1 ^ 3 - 

4 In this paper we use the Hill factor adopted by Henon & Petit 
(1986), which differs by 3 1//3 from the scaling used by some other 
authors. 
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advanced stages of dynamical evolution of the dispersion- 
dominated planetesimal population but it may fail in 
more general situations. 

It is thought (Kokubo & Ida 1998; Rafikov 2004; Gol- 
dreich et al. 2004) that at the very end of the oligarchic 
stage of planetary growth in the inner parts of the So- 
lar System, just before the transition to a chaotic final 
stage of planetary assembly (sometimes called the stage 
of giant impacts), planetesimal coagulation produced a 
number (several hundred) of protoplanetary cores with 
masses comparable to the mass of the Moon or Mars, i.e. 
~ (0.01 — 0.1) M®. These cores comprised a significant 
fraction of all the refractory mass of the disk and were 
well-separated in semi-major axis (typically by several 
Hill radii). 

The orbits of these cores are initially not expected to 
cross because their eccentricities are very small as a re- 
sult of efficient dynamical friction exerted on them by 
the residual population of small planetesimals. How- 
ever, with time the population of small planetesimals 
gets eroded by collisional grinding and accretion by cores 
and the strength of dynamical friction goes down. Dis- 
tant mutual gravitational perturbations between nearby 
cores then gradually increase their velocity dispersion, 
eventually allowing their orbits to cross, which leads to 
collisions between embryos and their growth into bigger 
bodies. 

Since initially the inclinations of the cores were almost 
zero distant perturbations cannot efficiently excite ver- 
tical motion of cores. Moreover, even though dynamical 
friction from the remaining planetesimals is no longer 
efficient in curbing the eccentricity growth of the cores 
it may still be strong enough to continue damping their 
inclinations. As a result, the protoplanetary cores are 
expected to reside in a very thin disk with i <C 1 all 
the way until the point when their orbits start to cross. 
When this happens one finds that i <C 1 < e, so that the 
condition i ~ e usually assumed in studies of planetesi- 
mal dynamical evolution is strongly violated. 

Thus, a vertically thin, dispersion-dominated planetes- 
imal disk can naturally arise in some circumstances. 
Since the collision rate of planetesimals is a sensitive 
function of their inclination - the smaller is inclination 
the higher is the collision probability and the faster is 
protoplanetary growth - it is important to know how 
much time the population of cores spends in the thin 
disk configuration after their orbits become crossing. If 
this time is sufficiently long then core masses could grow 
significantly by collisions even during the transient pe- 
riod when their inclinations have not yet increased. This 
possibility potentially may act to speed up the final as- 
sembly of terrestrial planets. 

A similar situation arises when one considers the tran- 
sition of the shear-dominated planetesimal population 
into the dispersion-dominated regime. Ida & Makino 
(1992) showed that a planetesimal population starting in 
the shear-dominated regime typically undergoes a phase 
in its dynamical evolution when i <C 1 < e (see their Fig. 
6 for illustration). This phase does not persist for very 
long but while it lasts the dynamics of the planetesimals 
are significantly different from the usually assumed case 
of i ~ e. 

These considerations give us a motivation to explore 



the dynamical regime i <C 1 < e in this work. The pa- 
per is organized as follows: in Sj2]we describe the equa- 
tions governing the velocity evolution of the planetesi- 
mals while in <j3] we analytically compute the scatter- 
ing coefficients entering these equations in the case of 
i <C 1 < e, and in 2] we compare our results with nu- 
merical calculations. In Sj5]we use our results to examine 
the velocity evolution of a population of protoplanetary 
cores with crossing orbits. In ^we provide comparison 
with other studies and discuss additional applications of 
our results. 

2. VELOCITY EVOLUTION. 

To understand the velocity evolution of planetesimals 
we consider two populations of planetesimals with masses 
mi and mi\ populations of different mass contribute lin- 
early to velocity evolution so it is sufficient to consider 
just two masses. We assume that for every planetesi- 
mal type ek,ik 1 and nik <C M*, k = 1,2, providing 
justification for using the Hill approximation. 

In this local approximation the Keplerian orbit of a k- 
th planetesimal type is described by the following equa- 
tions: 



Xk = h k -e k cos(t~T k ) 
3 

Vk 



Aft - ^h k t + 2esin(i - r fc ), (5) 



z k = i k sin(t - U! k ), 



(4) 
(5) 
(6) 



where x, y, and z are Cartesian coordinates in the lo- 
cal radial, azimuthal, and vertical directions centered at 
some reference stellocentric distance, h is the planetes- 
imal semi-major axis separation from the origin of this 
coordinate system, and A is a constant related to the 
origin of time t (measured in units of ll" 1 ). 
The relative motion of two non-interacting planetesi- 



mals in Hill units (f r 
by equations 



= (ri - r 2 )/aQui +fJ-2) 1/3 ) is given 



x r = h r — e r cos(i — r r ), (7) 
3 ~ 

y r = A r — -h r t + 2e r sin(t — r r ), (8) 

z r — i r sm(t — uj r ), (9) 

where e r , i r are the relative eccentricity and inclination 
of the planetesimals, h = (ai — 02)/a(/xi -I- /i2) 1 ^ 3 is the 
semimajor axes separation normalized in Hill units, and 
A r = (Ai — A 2 )/a(/i 1 + ^2) 1 ^ 3 - In the following we will 
drop the subscript "r" from all variables characterizing 
relative motion of planetesimals where it will not lead to 
confusion. 

Because of the mutual gravitational attraction relative 
orbital elements appearing in equations (0-© do not 
remain constant but change according to the following 
set of equations (Hasegawa & Nakazawa 1990; Tanaka & 
Ida 1996): 



dh 
~dt 

dX 
~dt 

de x 



dy'' 

d(f> d(j> 



dx 



dy 



. d<t> „ d<j) 
-suit— - 2cos£— , 
at ox ay 



(10) 

(11) 
(12) 
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de v dd> d<b 

— = cost— 2sinr— , 

at ox ay 

di x 

~~dt 

di», 



cost—, 
oz 



(13) 
(14) 
(15) 



~i = - sint 9i' 

where 

= -(£ 2 +y 2 + 5 2 )- 1 / 2 . (16) 

is the interaction potential. 

In this work we assume for simplicity that e and i of a 
fc-th planetesimal population have Gaussian distribution: 



ip(e k ,i k )de k di k 



de k di k 



■ exp 



2ai 



2aj 



A 7 ) 

where a Etk and are the dispersions of eccentricity 
and inclination of the fc-th population. Ida & Makino 
(1992) have found that a Gaussian distribution accu- 
rately describes the distribution of e and i found in 
direct N-body three-dimensional (3D) simulations of a 
large number of planetesimals gravitationally interacting 
in the dispersion-dominated regime. At the same time in 
their dispersion-dominated simulations of 2D disks Ida & 
Makino found more high-energy particles than a Gaus- 
sian distribution would predict. Despite this we still use 
distribution (|17p to represent velocities of planetesimals 
in thin disks as this is not going to strongly affect the 
velocity evolution but allows significant simplification. 

Our goal is to find how <7 e) i and cr^i of a population 
with mass mi varies in time as a result of gravitational 
interactions with planetesimals of mass 777,2 which have 
eccentricity and inclination dispersions <r ej 2 and 0^2 (for 
now we neglect other factors that may affect planetesi- 
mal velocities such as gas drag, inelastic collisions, and 
so on). General evolution equations for the case of dis- 
tribution (|17[) have been previously derived by a number 
of authors (Hornung et al. 1985; Ida 1990; Wetherill & 
Stewart 1993; Stewart & Ida 2000). Here we adopt a 
specific expression from Rafikov (2003): 



Ot 



= -nN 2 a 2 (m 

2 4 



/'2 



Mi + M2 



M2) 4/3 

//■2 



Mi 



M2 cr 2 1 



-H 2 



>e,2 



,(18) 



where N 2 is the surface number density of bodies with 
mass m2. 

Dimcnsionless stirring coefficients H\ 2 appearing in 
equation (fl"8)) are defined as 

H k = [ dedii> r {e,i)H k (e,i), k = l,2, (19) 



(20) 



H 1 = / dh\h\((AeY 



H 2 = / dh\h\((e- Ae)) T , 



(21) 



Here Ae is the change of e in the course of scattering, 
and (...)i- lU) = (4-7T 2 ) -1 J J drduj implies averaging over 



the relative orbital phases characterizing vectors e and i. 
Function t/v(e r , i r ) is the distribution function of relative 
e, i and can be shown (Stewart & Ida 2000; Rafikov 2003) 
to be given by 



ip r {e, i)dedi 



ede idi 



■ exp 



2a 



*2 



1 

25? 



(22) 



where ct 2 = (ct 2 1 + crf 2 )/(Ml + M2) 2/3 and af = (cr 2 ^ + 
(T ?2)/(/ i i+M2) 2 ^ 3 are the dispersions of relative eccentric- 
ity and inclination. Thus, functions H\ and H 2 represent 
scattering coefficients for a planetesimal population with 
a single value of relative eccentricity e and inclination i, 
while Hi and H 2 are these coefficients averaged over the 
distribution (|22| of e and i. 

The term inside the brackets in equation (|18j) propor- 
tional to Hi is called gravitational stirring (Rafikov 2003) 
while the term proportional to H 2 is called gravitational 
friction and is different from the "dynamical friction" 
used by other authors (Stewart & Ida 2000; Ohtsuki et 
al. 2002). 

Equation (fP8|) also describes the self-stirring of pop- 
ulation with mass mi if one changes the subscript "2" 
to "1" in its right hand side (which makes expression in 
brackets equal to (Hi + 2H 2 )/A). For a continuous dis- 
tribution of planetesimal masses equation (|18p should be 
generalized by integrating the right hand side over the 
mass spectrum. 

Equations analogous to (|18|) can be written for the in- 
clination evolution by replacing "e" by "i" everywhere in 
equations (fT8| - ((2T|) and using scattering coefficients Ki. 2 
and K\ t 2 defined analogously to expressions (fT9|) - ((2T|) in- 
stead of Hi y2 and #1,2)- 

3. SCATTERING COEFFICIENTS. 

System ([T8|) - ([2T]) is a rather general set of equations 
derived for a Gaussian distribution of orbital elements 
(|17|) . It allows one to determine how o~ e>k and ai yk (k = 

1,2) evolve in time once coefficients Hi 2 and Kip are 
known as functions of a e and (7j. These coefficients have 
been previously calculated by a number of authors in the 
dispersion-dominated case under the assumption that i ~ 
e. However, as we demonstrate shortly, these calculations 
become invalid once i gets sufficiently small while e> 1. 
Thus our next step is to rederive scattering coefficients 
in the case of i -C e from the first principles. 

In doing so we will adopt an approach previously de- 
veloped by Ida et al. (1993) for the dispersion-dominated 
regime. According to this method (1) the approach tra- 
jectory of one planetesimal far from another is repre- 
sented as a straight line and (2) the effect of perturba- 
tions from the central mass on the gravitational scat- 
tering of the two planetesimals is neglected. There is 
also an implicit assumption that (3) the scattering coef- 
ficients are dominated (or at least strongly contributed 
to) by those planetesimals whose trajectories pass very 
close to the perturber. If the latter assumption is ful- 
filled then the other two are quite natural. Indeed, the 
first approximation works well because most of the per- 
turbation to the orbit of the passing particle occurs near 
the point of closest approach of the interacting bodies. In 
this region the curvature of the planetesimal trajectory 
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Dispersion-dominated regime 



e ^ 



Thin disk regime 



Fig. 1. — Schematic illustration of different regions in the e- 
i phase space, showing the shear- and dispersion-dominated re- 
gions and the thin-disk dynamical regime (shaded), which is also 
dispersion-dominated (e > 1). 



caused by epicyclic motion can be neglected, justifying 
the straight-line simplification. The second approxima- 
tion works because the most significant contribution to 
scattering is due to trajectories passing very close to the 
perturber, within its Hill radius, where the influence of 
the third, central body can be disregarded. 

When e ~ i ;» 1 the flux of approaching planetes- 
imals in the vicinity of any given perturber is roughly 
uniform on scales ~ Rh, and every decade in the ini- 
tial impact parameter of interacting bodies contributes 
roughly equally to the scattering coefficients (Binney & 
Tremaine 1987). This gives rise to appearance of the 
so-called Coulomb logarithm, In A, in expressions for the 
scattering coefficients. The argument A is roughly the 
ratio of the maximum impact parameter Z max ~ iRn, 
beyond which the density of approaching planetesimals 
is non-uniform, to the impact parameter 

Rh (23) 



I, 



e 2 + i 2 ' 



at which the trajectories of incoming planetesimals ex- 
perience large-angle deflection. 

Previous calculations of the scattering coefficients in 
the dispersion-dominated regime assumed that lmax 3> 
l m in meaning that In A > 1. In this case the scattering 
calculation for closely approaching orbits with impact pa- 
rameters ~ lmin Rh which can be done analytically 
approximates quite well (with logarithmic accuracy) the 
full scattering coefficients so that the aforementioned as- 
sumption of the dominance of close encounters for scat- 
tering coefficient calculation is roughly fulfilled. Clearly 
lmax 3> lmin requires that i 1/e 2 , which is essentially a 
condition for the standard expressions for the dispersion- 
dominated scattering coefficients to be valid. 

In this work we look at the opposite extreme, namely 
a situation when 



< 



i crit = e 2 < 1, e > 1, (24) 

(see Figure Q] for illustration) . When this condition is 
satisfied the assumption of a uniform distribution of ap- 
proaching planetesimals around the scatterer does not 
hold even for I ~ l m in- Then a new calculation of scat- 
tering coefficients is needed. 

In Appendix A we present the details of such a calcu- 
lation which makes the following set of assumptions: (1) 



the planetesimals move at high relative velocities which 
allows us to use a two-body scattering approximation, 
(2) the planetesimal velocities change only during close 
approaches, which are possible only when h < e, (3) the 
gravitational interaction between planetesimals at large 
separations is neglected, and (4) after changing as a re- 
sult of the encounter with the scatterer the planetesimal's 
orbital elements do not change further. These simplifi- 
cations allow us to derive the following set of expressions 
for the integrands of the scattering coefficients, see equa- 
tions P0 ]) -([2T ]) : 



((Ae) 2 )_ 
(e • Ae) WjT 
<(Ai) 2 ).,. = 
(i • Ai)o),T = 



20 



\h\Ve?-h 2 ' 



4 e 2 



3 v\h\Ve 2 — h 2 
v 5 i 2 



2 ,~572 



\h\Ve 2 -h 2 



3 r. 



v\h\Ve 2 -h 2 



(25) 
(26) 
(27) 
(28) 



Note that these expressions do not contain a Coulomb 
logarithm and do not suffer from the ambiguity related 
to the choice of minimum and maximum values of im- 
pact parameter l m in and l m ax typical for the 3D case. 
The physical reason for this lies in the fact that in the 
limit (f2"4")) the scattering coefficients are dominated by 
trajectories with impact parameters I ~ l m im i-e. those 
leading to large-angle scattering. Thus, integrals over dl 
appearing in the calculation of the scattering coefficients 
are mostly contributed to by I ~ l m m in the quasi-2D 
case 5 , so that values of I much larger and much smaller 
than l m in affect the coefficients only weakly. This is very 
different from the 3D case, in which trajectories expe- 
riencing weak scattering provide significant contribution 
to the scattering coefficients. 

Integrating these equations over \h\dh from to e (lim- 
its within which a given planetesimal can experience a 
close encounter with the scatterer) and substituting into 
equations (f2"0")) - (|2ip we arrive at the following expressions 
for the scattering coefficients corresponding to fixed e and 
i and averaged over the phase angles r and u>: 



Hi 



H 2 = 



16.147 e r , 




0.864 i 2 e 5 , 



Kn 




-2.875 ^. 



(29) 
(30) 

(31) 
(32) 



Finally, averaging coefficients (|2"9"|) - ([3"2")) over the Gaus- 

5 Note that this property makes the assumptions adopted in our 
calculation quite robust. 
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sian distribution (fT7|) one finds that 

H x = E *e,r « 20.237 a e , r , (33) 



7.207 CT e;r (34) 




K 2 = — 5L_K -^ U^w -7.207 ^ (36) 

These expressions represent the behavior of scattering 
coefficients in the limit (1241). 



4. COMPARISON WITH NUMERICAL RESULTS. 

To check our analytical results we ran a series of numer- 
ical calculations. The latter are the Monte-Carlo com- 
putations of integrals in equations (fl9l) - pT|) with Ae r , 
Ai r obtained by direct integration of equations (fTT)|) - (fl"5|) . 
Equations for the evolution of orbital elements have been 
integrated using fourth-order Runge-Kutta method with 
adaptive step size control (Press et al. 1992). Conserva- 
tion of Jacobi constant has been routinely monitored and 
this integral of motion has been found to vary during the 
calculation by at most one part in 10 5 for a very small 
number of orbits. In the majority of calculations the Ja- 
cobi constant has been conserved to relative accuracy of 
better than 10~ 10 . 

The orbits used in the Monte-Carlo evaluation of inte- 
grals have been drawn from the distribution of initial or- 
bital parameters appropriate for each particular scatter- 
ing coefficient. When computing (e- Ae) W)T , ((Ae) 2 ) WjT , 
etc. we select a set of values of e, i, and h, and draw r 
and w randomly from a uniform distribution between 
and 2ir. When computing H12, ^1,2 we also draw h from 
a uniform distribution between —Lh and Lh, while keep- 
ing e, i fixed. Finally, to compute Hi 2, -Ki,2 we draw 
e x , e y , ixi iy randomly from a Gaussian distribution (I17| 
with given dispersions d e and Oi, while h is drawn from 
a uniform distribution between —Lh and Lh- 

In all our of calculations of -Hi, 2, Ki,2 we use 

L h = 10 + 4c, (37) 

to ensure that even orbits with h > e are properly ac- 
counted for. When computing -Hi. 2, If 1,2 we use the same 
prescription for Lh but with a e replacing e. We adopted 
the following prescription for the number of orbits used 
for evaluating scattering coefficients: 

iV = 5x 10 4 [10+ (1 +i)(l + e) 2 ] . (38) 

Thus, to compute scattering coefficients for e = 15 we 
run around 13 million scattering calculations. For H\ 2, 
K\ t 2 we used the same prescription with a e , di replacing 
e, i. 

In Figure [5] we present the results of calculation of 
Hi, 2, Ki^/i 2 as a function of e, for several values of 
i = 10 _1 , 10~ 2 , 10~ 4 , 10~ 6 , together with our analytical 
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Fig. 2. — Results of numerical calculation of scattering coeffi- 
cients (a) Hi, ( b) H 2 , (c ) Ki/i 2 , (d) K'2/i 2 , compared with analyt- 
ical predictions 129M321 . represented by solid lines. Values of coef- 
ficients are shown as functions of e for i = 10~ 1 , 10~ 2 , 10 — 4 , 10 — 6 
(see legend in panel (b) for associating different dot styles with 
particular i). 
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predictions (|29)) - (f32]) shown by solid lines. Values of K1.2 
are scaled by i 2 to simplify comparison of curves corre- 
sponding to different i. 

Figure [^i demonstrates rather good agreement be- 
tween analytical and numerical results for the gravita- 
tional stirring coefficient Hi almost everywhere in the 
considered range. Agreement at very small values of e is 
likely a coincidence since at e sa 1 a transition to a shear- 
dominated scattering should occur which invalidates our 
assumption of e S> 1 . At higher values of e curves corre- 
sponding to different i generally line up with the analyt- 
ical results quite well except for the noticeable deviation 
of i — 0.1 curve from analytical result (|29|) which starts 
around e « 3 and becomes stronger as e grows. This de- 
viation is expected since our analytical results are valid 
only for i satisfying the constraint (|24[) . For i = 0.1 this 
means that agreement with the analytical result is ex- 
pected only for e < i~ x l 2 « 3, in good correspondence 
with Figure [2l This point is reinforced by observing the 
deviation of i = 10~ 2 results from the analytical curve 
that starts at e w 1 w (10 -2 ) 1 / 2 , i.e. again agrees with 
the constraint (|2g|). Points for i = 10 -4 and 10~ 6 fall 
on top of each other in the whole range of calculation 
as they should. They lie somewhat below the analytical 
prediction for large values of e, which we do not have a 
good explanation for. 

The same applies very well to the results for the grav- 
itational friction coefficient H2 shown in Figure^. The 
only slight difference is that the influence of the shear- 
dominated regime is more pronounced for this scattering 
coefficient, as H2 settles onto the analytical result (|30p 
only at e ~ 2.5. Thus, the scattering coefficients based on 
eccentricity changes agree with theory quite well within 
the range of applicability of the analytical results. 

We now turn to coefficients Kx,2 which are based on 
changes in inclination. As one can see from Figure[2p, the 
shear-dominated regime affects stirring coefficient K\ for 
e < 2.5. As expected from our previous discussion, K\ 
strongly deviates from the analytical prediction starting 
at around e ~ 3 for i = 0.1 and at around e ~ 10 for 
i — 10~ 2 . However, the results for i = 10~ 4 and i — 10~ 6 
generally do not fall on top of each other as one would 
expect given that the quadratic scaling of K% t 2 with i 
has been removed in Figures [2k, d. Moreover, the values 
of K\ clearly deviate quite strongly from the analytical 
prediction (|31|) . sometimes by three orders of magnitude, 
without any recognizable regular pattern. 

The numerical results for the gravitational friction co- 
efficient K2 as compared with theory are even more sur- 
prising, as Figure[2ji demonstrates. Here, for small values 
of e < 5, K2 systematically increases with e, while ana- 
lytical result (|32[) predicts that K% should be a decreasing 
function of e. At larger e numerically computed values of 
K2 exhibit significant scatter in a chaotic fashion. To be 
fair, one should note that some of the theoretical expec- 
tations are confirmed by numerics even for K2: a curve 
for i = 0.1 again diverges from other curves correspond- 
ing to smaller values of i at e ~ 3. 

In Figure [3] we look at the behavior of scattering coef- 
ficients as functions of i for a fixed value of e. The main 
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goal of these plots is to illustrate the transition between 
the thin and thick disk regimes of planctcsimal scattering 
occuring at i crit ~ e~ 2 . The rather good accuracy of our 
analytical results can be clearly seen in the behavior of 
Hi t 2 and even K\ (the situation is less clear in the case of 
K2Y the behavior of the scattering coefficients changes 
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dramatically at i w 10 2 for e = 10 and at i w 0.05 for 
e = 4. The long-dashed lines in Figure [3] illustrate the 
scaling of the scattering coefficients with i in the thick- 
disk regime (Stewart & Ida 2000), and show good agree- 
ment with our numerical results when the condition (|24| 
is violated. It is clear from Figures [5^,b that our thin- 
disk theory describes the behavior of eccentricity-based 
coefficients quite accurately even for e = 4, which 
is not very far from the shear-dominated regime. 

However, from Figures [3j:,d one sees once again that 
the inclination-based coefficients Ki t 2 deviate from ana- 
lytical predictions. Already for e = 4 coefficient K\ ex- 
hibits stochastic variations as a function of i by a factor 
of order unity. At e = 10 these variations become quite 
dramatic and exhibit an increasing trend with decreasing 
i. This is rather surprising since one expects analytical 
theory to work better for very small values of i, when the 
condition (pM)) is satisfied by a large margin. This clearly 
indicates that the theory is missing some important in- 
gredient, a conclusion which is additionally reinforced by 
Figure [3ji demonstrating rather poor agreement between 
numerical and analytical values of K^. 

In Figures BJ [5] we show the behavior of scattering co- 
efficients Hi,2) -^4,2 averaged over the Gaussian distribu- 
tion of e and i. As expected, all the major features of 
Hi,2, K iy 2 discussed above are preserved in these plots, 
although the overall agreement with theory is addition- 
ally spoiled by the fact that numerically computed H\2, 
K1.2 represent a Gaussian convolution of -ffi^, ^1,2 over 
an extended range in e and i, and not everywhere inside 
this range are the basic assumptions (e.g. dispersion- 
dominated scattering) of our analytical theory fulfilled. 
In particular, numerical coefficients are affected to some 
extent by shear-dominated scattering events, not ac- 
counted for in our theory. Also, at high di ~ 0.1 — 10~ 2 
a significant fraction of numerically integrated scatter- 
ing events had values of i ~ e corresponding to thick 
disk scattering, for which the behavior of coefficients is 
different from our theory; see Figured) 

Based on the results presented in Figures [2][5] we con- 
clude that analytical theory explains quite well the be- 
havior of scattering coefficients based on changes of e, 
while it provides a rather poor fit to the numerically 
determined behavior of the inclination-based scattering 
coefficients. There may be several reasons for this dis- 
crepancy, some of which are listed below. 

1. It may be that the discrepancy arises when we in- 
tegrate the phase-averaged coefficients (e • Ae) WjT , 



etc. over h to obtain Hi 2, etc.; see 



definitions (|2H)l and (|2~lj) . In particular, encounters 
with h > e neglected in our analytical work may 
provide an important contribution to the numeri- 
cally computed rates. 

2. The two-body approximation used in our analytical 
calculations does not work well. 

3. Our assumption of a single close scattering per ap- 
proaching orbit may be faulty, as the scattered 
planctcsimals may have orbital parameters allow- 
ing them to experience additional close approaches 
with the scatterer. 
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Fig. 6. — Plots of phase-averaged scattering coefficients (a) 
{(Ae)\,r, (b) (e • Ae) u , T , (c) {(Ai) 2 ),, T , (d) (i • Ai)^ T , as 
functions of e for a fixed value of h = 10 and several values of 
i = 0.1, 10~ 2 , 10~ 4 , 10 -6 ; see legend in panel (b). Note the rapid 
decay of scattering coefficients for e < h. Dotted lines show ana- 
lytical predictions for e > h. 



4. Our theory assumes that the changes in planetesi- 
mal orbital elements occur only during the close ap- 
proach, when the planetesimal separation is < Rh , 
while in reality it may be that the more distant 
interactions between planetesimals at separations 
> Rh also play an important role. 

We devote the rest of this section to exploring these 
possibilities. 

4.1. Integration over h. 

To figure out whether the aforementioned discrepancy 
between the analytical and numerical inclination-based 
scattering coefficients can be caused by the integration 
of the phase-averaged coefficients over h we look at the 
behavior of the phase-averaged coefficients. In Figure 
[S] we present their scaling with e for a fixed value of 
ft, = 10 and several values of i. Similarly, in Figure [7] 
these coefficients are shown as functions of h for fixed e 
and the same values of i. Based on these plots we can 
make several conclusions. 

First, when e > h analytical predictions for (e • Ae) UiT 
and ((Ae) 2 ) W T fit our numerical results quite well. Co- 
efficients computed for i = 0.1 deviate from theory be- 
cause, as previously described, they do not correspond 
to the thin disk scattering regime. At the same time 
(i • Ai)^ iT and ((Ai) 2 ) w r are still significantly different 
from theory and exhibit rather erratic behavior. 

Second, the values of all scattering coefficients corre- 
sponding to e < h are much smaller than their values for 
e > h. In the latter case planetesimals can experience a 
close approach, while in the former this is not possible, 
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Fig. 7. — Same as Figure [7] but with phase averaged coefficients 
plotted as functions of h for e = 10. 



and changes of orbital elements are much weaker than in 
the latter case. As a result, contribution of orbits with 
e < ft to scattering coefficients is very small. In fact, one 
can see from Figure [7] that for e = 10, ((Ae) 2 ) w r com- 
puted at ft = 10 (close encounters possible) and h = 11 
(close encounters not possible) differ by more than 2 or- 
ders of magnitude. The same is true for (e ■ Ae) WiT , 
while for ((Ai) 2 ) w , T and (i • Ai) w , T this difference is 6 
and 3 orders of magnitude respectively. This very well 
illustrates our point made in that trajectories expe- 
riencing large-angle scattering (possible only for e > ft.) 
strongly dominate scattering coefficients in the case of a 
thin planctesimal disk. 

To summarize, the results displayed in Figures [SH 
make it clear that the well-separated orbits with e < ft, 
for which strong scattering is impossible, do not con- 
tribute much to the scattering coefficients. All the dis- 
crepancy between the theoretical and numerical values 
of K\ t 2 and K\ 2 is already present in the corresponding 
phase-averaged coefficients, and is not introduced by the 
integration of the phase-averaged coefficients over only a 
finite range of ft, \h\ < e. 

4.2. Accuracy of the two-body approximation. 

Next we consider whether the two two-body approxi- 
mation adopted in our analytical calculations is valid for 
the case of thin disk scattering. 

Previously, Tanaka & Ida (1996) have compared nu- 
merically computed changes of orbital parameters re- 
sulting from gravitational scattering with analytical pre- 
dictions derived in the two-body approximation. They 
found good agreement between the two, except for the 
narrow regions of the initial epicyclic phases in which 
orbital parameters evolved in a chaotic manner, if (a) 
ft > 2, (b) the encounter velocity vq > 4 and (c) a small 



shift in the initial epicyclic phases r and lo is introduced 
to match analytical predictions. Tanaka & Ida (1996) 
have compared only Aft calculated by both methods for 
i = 0, and also the changes of other orbital elements 
for i ~ e > 1. None of these cases corresponds to the 
regime of thin disk scattering considered in this work al- 
though the former does describe quite well the variation 
of the eccentricity-based scattering coefficients. For that 
reason we ran our own calculations with initial orbital 
parameters selected to correspond to the thin disk case. 

In general we find good agreement with the conclu- 
sions of Tanaka & Ida (1996), as shown in particular in 
Figure [8] where we display the changes of various orbital 
elements resulting from gravitational scattering as well 
as the minimum approach distance between the scatter- 
ing bodies l m i„. In making this Figure we have slightly 
shifted analytical curves (shown as dotted lines) in r by 
At = —0.05 to make them better match numerical re- 
sults (shown as solid curves). In practice such a shift 
of the orbital phase arises due to the distant interaction 
between the planetesimals as they approach each other 
(Tanaka & Ida 1996). Only one interval of r in which 
strong scattering is possible is shown, 0.3 < r < 0.39; 
another one exists at 5.97 < r < 6.05, in accordance 
with the discussion in Appendix, where the existence of 
two values of r for which close encounters are possible 
for e > ft is stated. 

One can deduce from Figure [8] that analytical curves 
follow the numerical results quite well for the majority of 
values of r except for the two narrow ranges of r, namely 
0.32 < r < 0.33 and 0.35 < r < 0.355. Inside these inter- 
vals orbital elements experience strong chaotic variations 
as t changes, with Aft, Ae^, Ae y deviating from analyt- 
ical prediction by a factor of order unity, while Ai x , Ai y 
differ from theory by several orders of magnitude (off 
scale on these plots) ! Although these deviations are very 
significant we will show next that they are not caused by 
the failure of the two-body approximation. Thus, use of 
the two-body approximation cannot explain the discrep- 
ancy between the numerical and analytical inclination- 
based scattering coefficients. 

4.3. Single-scattering approximation. 

Our calculation has always assumed that changes of or- 
bital elements arising during a scattering event are final. 
In reality there may be a situation when a post-scattering 
orbital elements are such that they cause another close 
approach between planetesimals, leading to additional 
variation of orbital elements. And this may happen not 
just once for a given incoming orbit. Such multiple scat- 
tering events are very typical for planetesimals scattering 
in the shear-dominated regime but their importance in 
the high- velocity case is not very obvious. 

To see that multiple scattering is indeed possible even 
in the dispersion-dominated regime we take a closer look 
at Figure[8l where we plot Aft as a function of r. One can 
see that chaotic orbits strongly deviating from analytical 
prediction (shown as a dotted line) exist almost solely 
in those regions where Aft predicted by theory happens 
to be f» —ft (Aft w —5 w —ft in the case displayed in 
this Figure). This is not a coincidence, and what re- 
ally happens is the following. First, planetesimals scat- 
ter and their orbital elements change in full agreement 
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Fig. 8. — Changes of relative orbital parameters (a) (Ai) 2 , 
(b) i ■ Ai, (c) (Ae) 2 , (d) e • AS, (e) Ai*, (f) Ai y , (g) Ae x , (h) 
Ae y , (i) Ah, and (j) the minimum separation l m i n of planetesimals 
plotted as a function of their relative horizontal epicyclic phase r 
(at large separation prior to scattering). Initial orbital parameters 
corresponding to this calculation are shown at the top of the plot. 
Only an interval 0.3 < r < 0.39 in which strong scattering takes 
place is displayed. Solid curves show the numerical results while 
the dotted lines are the analytical predictions. Note the chaotic 
variation of orbital parameters for 0.32 < r < 0.33 and 0.35 < r < 
0.355. 



with analytical theory. This means, however, that post- 
scattering h ~ 1 and the guiding center of the planetesi- 
mal orbit is now moving very slowly with respect to the 
scatterer, while eccentricity is still quite high. Right af- 
ter scattering planetesimals are very close to each other, 
and Keplerian shear does not allow their guiding centers 
to recede very far because h is small, so that after one 
orbital period planetesimals may closely approach each 
other again and experience another scattering. This sec- 
ond scattering may or may not dislodge them from close 
proximity of each other but it will certainly affect their 
final orbital elements, explaining the deviation of Ah and 
other orbital elements from theoretical predictions when 
Ah « —h. Thus, we conclude that 

• Chaotic variations of orbital elements result from 
multiple scatterings of planetesimals in the course 
of close encounter rather than from the failure of 
the two-body approximation. 

• Phase intervals where theoretical Ah « —h are nat- 
urally occupied by chaotic orbits (there are also 
other possibilities for producing multiple scatter- 
ing orbits, see below). 

These observations greatly help in explaining the puz- 
zling results for the inclination-based scattering coeffi- 
cients. From pure geometry it is clear that the highest 
inclination i\ which a high- velocity particle with initial 
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(39) 

This is easily seen in Figures |S^,f which show that max- 
imum ix ~ 10 -2 for io ~ 10 -4 : since vq — (e 2 — 
(3/4)/? 2 ) 1 / 2 w 4.2 one should expect maximum i\ ~ 
10~ 4 x 4.2 3 « 0.007, very close to what we find in re- 
ality outside of the region of chaotic orbits. 

At the same time Figures |5j3,f show that chaotic orbits 
often exhibit final i much larger than predicted by equa- 
tion (f39|) . This, of course, is naturally explained by the 
fact that chaotic orbits result from multiple scattering. 
Every scattering of a high-velocity orbit can potentially 
increase inclination by a factor Vq ^> 1 (the approach 
velocity of planetesimals does not change very strongly 
after multiple scatterings and is still ~ vq] this can be 
understood from the conservation of Jacobi constant). 
Thus, after n scatterings maximum possible inclination 
would have been i n ~ ioVQ U , except that in practice i 
cannot exceed Vq. The highest final i that we could find 
for the parameters of Figure [8] is w 0.7 but one has to 
keep in mind that orbits in the chaotic region exhibit 
quasi-fractal behavior in that the denser is the grid in 
r used for computing Ai the richer the behavior found. 
Thus, we could have easily missed orbits with even higher 
final i. On theoretical grounds we expect the maximum 
i in this Figure at the level of 10~ 4 x 4.2 3x2 m 0.5 if only 
n = 2 close scatterings have taken place. 

In Figure [9] we illustrate a multiple scattering event 
for an orbit with initial parameters e = 15, i = 10 -6 and 
h w 12.6 (in this case vq m 10.2). From Figure [9Jd,c one 
can see that as a result of the first scattering h becomes 
very small while i jumps up by - 10 2 . After that the 
planetesimal loops around its scatterer for several orbital 
periods, as illustrated in Figures [9ji-f, until the second 
strong scattering takes place, resulting in final h m 20. 
This allows the planetesimal to leave the vicinity of its 
scatterer. During the second scattering i is again boosted 
up by more than two orders of magnitude, resulting in 
final i w 0.03. This is much larger than 10~ 6 x 10.2 3 » 
10 -3 — the maximum i\ one would expect from single 
scattering. 

It now becomes much easier to understand the erratic 
behavior of scattering coefficients K\ <x in Figure [2] In 
particular, multiple scattering orbits very strongly affect 
K\ since Az enters the calculation of this stirring coef- 
ficient in a second power. Because of that, even though 
chaotic orbits arise for only a small subset of horizontal 
epicyclic phases they affect the value of numerically de- 
termined K\ very strongly. According to equation (127|) 
((Ai) 2 ) WjT w 6.5 x 10~ 10 for the values of e, h and i used 
in making Figure [H At the same time a single orbit as 
displayed in Figure [9] has (Az) 2 ~ 10~ 4 , more than 10 6 
times higher than the theoretical phase average of this 

6 Highest inclination results from scattering by fa 7r/2, which 
requires impact parameter of incoming trajectory to be I ~ , 
at initial vertical separation of order i^Ru- The final velocity 
of the receding planetesimal is Co and from simple geometry its 
vertical component (which is equivalent to inclination in Hill units) 
is v X (io/l) ~ ioVg. 
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Fig. 9. — (a-c) Variation of the relative orbital elements of 
two planetesimals in the course of a multiple scattering event, for 
initial orbital elements indicated at the top of the plot. Evolution 
of (a) the relative distance between the bodies r, (b) their relative 
inclination i, and (c) their e (solid line) and h (dotted line) are 
shown. This trajectory exhibits two strong scattering events in 
each of which i gets boosted up by two orders of magnitude, (d-f) 
Trajectory of relative motion in the course of scattering shown (d) 
in the y — x coordinates and (e) in y — a coordinates. A zoomed 
in version of panel (e) is shown in panel (f) to better illustrate 
the complexity of vertical motion before the final scattering causes 
planetesimals to recede from each other. In panel (d) we show both 
the instantaneous position of planetesimal that is being scattered 
(solid line) and the trajectory of its guiding center (dotted line). 
See text for more details. 



quantity. It is thus not surprising that even though we 
have used a very large number of orbits in computing 
scattering coefficients (according to the prescription ([38]) 
our calculation of K\ for e = 15 used 13.3 million orbits) 
chaotic orbits still affect them quite significantly. 

We can now explain why the scatter in numerical val- 
ues of K\ and the deviation from analytical prediction 
both become stronger as i decreases: the maximum pos- 
sible value of i resulting from scattering is always limited 
from above by i ~ Vo, so that the maximum stochastic 
(Ai) 2 ~ Vq, independent of initial i. However, the ana- 
lytical value of K\ cx i 2 , so that the ratio of analytical 
K\ to the numerical one increases as i decreases. 

It is not even clear that our calculation of K\ in Figure 
[5] has converged — one cannot guarantee that increasing 
the number of orbits would not increase even more the 
number of extremely chaotic orbits with very large Ai, 
which would then dominate the calculation. The only 
thing that argues against this scenario is the saturation 
of final i at the level of vq even for very large number 
of repeated scatterings. Nevertheless, until we under- 
stand how much of the phase space volume corresponds 
to chaotic orbits with very large Ai we cannot draw a fi- 
nal conclusion about the convergence of K\ and we leave 
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Fig. 10. — Same as Figure [9] but for a different choice of orbital 
elements indicated at the top of the plot characterized by small 
initial h. Note an exponential growth of i by several orders of 
magnitude in panel (b). In panel (d) we plot only a small fraction 
of data during the scattering event as dots to better illustrate the 
underlying structure. See text for more details. 



this subject for future investigation. Paradoxically, the 
agreement between the numerical and analytical results 
may be better if one uses smaller number of orbits in 
numerical calculation of since then the chance of 
randomly picking one of the high-Ai, multiple scattering 
orbits is also smaller. 

In the course of our investigation we have also found 
that orbits with e> 1 and /[ > 1 (like the one shown 
in Figure [9]) are not the biggest contributors to chaos 

in Ki 7 2- It turns out that inclination-based scattering 
coefficients arc most strongly affected by orbits with 
e ^> 1 and h ~ 1, i.e. orbits which are initially close 
to the separatrix between the horseshoe and passing or- 
bits. An example of planetesimal scattering correspond- 
ing to this case is shown in Figure [TOl for initial h w 1.02, 
e = 10, i = 10~ 6 . This event is characterized by a very 
long time interval, more than 10 3 orbital periods, during 
which planetesimals stay close to each other. They es- 
sentially form a temporary distant satellite system (note 
that the distance between planetesimals is larger than 
Rh), which slowly evolves in time. Figure [9)d demon- 
strates that during this temporary capture h oscillates 
around zero not allowing planetesimals ro recede from 
each other. Their relative inclination increases exponen- 
tially (with rather long time constant) by 4 orders of 
magnitude in an orderly fashion. Finally a strong scat- 
tering event occurs, which boosts up i by ~ 10 2 and 
dislodges the planetesimal from its scatterer's vicinity. 
For this event (Ai) 2 ~ 1, while a single scattering cal- 
culation would predict ((Ai) 2 ) w , r w 6.7 x 10~ 9 for the 
values of e, ft, and i shown on top of Figure [11] As a 
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2x10-" H( a ) 



^10- 



Fig. 11. — Same as Figure llOl but for a different choice of orbital 
elements indicated at the top of the plot. Note that while after 
the first scattering h is not particularly close to unity, multiple 
scattering still takes place. See text for more details. 



result, this single scattering event completely determines 
the calculation of K\ . 

It is worth pointing out here that multiple scattering 
in general does not require Aft. « — h and Figure [11] il- 
lustrates this statement. This Figure shows a double 
scattering event for initial h = 5, e = 12, and i = 1CP 4 . 
As can be seen in Figure [Tib after the first strong scat- 
tering event h ps 3.5 and Keplerian shear ensures that 
the bodies will not stay close to each other for very long. 
However, before the planetesimal leaves its scatterer's 
vicinity its epicyclic motion brings it back into their mu- 
tual Hill sphere where another scattering event occurs. 
This type of multiple scattering event does not affect co- 
efficients if 1,2 nearly as much as events with A/i« —h. 

To summarize, multiple scattering explains the dis- 
crepancy between the analytical and numerical results 
for K\ and stochastic scatter at high e in values of K% 
quite well. However, this explanation does not work so 
well for the systematic deviation of the numerical K.% 
from the analytical one clearly seen in Figure [2ji for vir- 
tually all values of e: even at small e, when the stochastic 
scatter is small, K% increases contrary to theory. Note 
that while the lower envelope of K\ for a given i agrees 
quite well with the analytical prediction (and multiple 
scattering explains the remaining stochasticity), this is 
clearly not true for K^. 

4.4. Distant interaction. 
To understand the systematic deviation of the numer- 
ical Ki from the analytical prediction ([32]) we first note 
that coefficient (i • Ai) W T used in calculating of Ki also 
systematically differs from the analytical prediction given 
by equation ([28]k see Figures [6ji, [Tji- This may seem 
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Fig. 12.— (a-d) Variation of (a) (Ai) 2 , (b) i ■ Ai, (c) Ai y , and 
(d) Ai x with vertical epicyclic phase ut for a particular set of ini- 
tial orbital parameters (shown on top of left column) . Theoretical 
predictions are shown with the dotted line while numerical results 
are in solid, (e-g) Plots of (e) (Ai) 2 and (f) i ■ Ai averaged over 
u) as functions of r for a particular set of initial orbital parameters 
(shown on top of right column). Solid lines arc numerical results, 
dotted lines are analytical predictions. Panel (g) shows the ratio 
of the numerically computed average of i ■ Ai over u) to the ana- 
lytical prediction for the same quantity. Panels (f) and (g) clearly 
demonstrate that theory underpredicts (i ■ Ai)^ by about an order 
of magnitude. 



surprising since in Figure [SJd the numerically determined 
i • Ai follows quite closely analytical prediction as a func- 
tion of horizontal phase r. Chaotic variations of i • Ai 
due to multiple scattering noticeable in this plot within 
narrow intervals of t, cannot explain the systematic dis- 
crepancy for K% - they can only be responsible for the 
random scatter in the calculation of K.2- However, Fig- 
ure [5}} does not show how i- Ai depends on oj (this Figure 
was made for a single value of u>) and this dependence 
turns out to be very important. 

In Figure [T2"R-d we display the dependence of (Ai) 2 , 
i ■ Ai, Aij,, and Ai x on lu for a fixed r. The value of 
r = 0.347 is chosen in order to avoid intervals of chaotic 
variation of the orbital elements 7 in order to isolate the 
subsequent analysis from the effects of multiple scatter- 
ing. The general agreement of the numerical and ana- 
lytical curves, including those of i • Ai, is quite good in 
this Figure. However, to calculate (i • Ai) W:T we need to 
integrate i ■ Ai over uj and it is quite obvious from Figure 
WBp that i • Ai is very close to a pure sinusoid. Its oj- 
average (i • Ai) w should then be very small and strongly 
dependent on the deviations of i • Ai from a pure sinu- 
soid. As a result, if the numerical and analytical values 

7 See Figure [8] whi ch i s made for the same set of initial orbital 
parameters as Figure [12] 
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of i ■ Ai deviate from the sinusoid differently, one can get 
a significant discrepancy between theory and numerical 
calculation. 

This is exactly what is going on as we demonstrate 
in Figure [T2b.f where we plot ((Ai) 2 ) w and (i • Ai) w as 
functions of r. One can see that averaging over ui does 
not affect the agreement between numerical and analyti- 
cal (Ai) 2 seen in Figure f!2k because it is an intrinsically 
positive quantity. But w-averaging of % ■ Ai does lead to 
a dramatic difference between analytical and numerical 
(i • Ai) u in Figure [127 for the reason we just described. In 
Figure WZk we show the ratio of (i • Ai) w determined by 
the two methods, and one can clearly see that the numer- 
ical result significantly exceeds the analytical one (both 
in regions of chaotic and orderly behavior of orbital pa- 
rameters) , in agreement with the fact that the numerical 
K% is systematically higher than the analytical K2, see 
Figure Eli. 

What may produce such a difference in scaling of an- 
alytical and numerical i ■ Ai with u)l In Appendix [Bl we 
provide a simple calculation of (i • Ai) w allowing for a 
small but non-zero difference Sujdist of the relative verti- 
cal epicyclic phase of interacting planetesimals u> at large 
separation and right before the close encounter. Such a 
phase difference in u> arises because of the distant interac- 
tion between planetesimals prior to their encounter and 
is analogous to the shift in horizontal phase r which was 
invoked in Figure [5] to better match analytical and nu- 
merical results (see also Tanaka & Ida 1996). We show in 
AppendixlBlthat although |<5oJdist | is expected to be small 
its effect on the calculation of (i • Ai) w is very important 
(and dominates this calculation) as long as Uol^distl > 1. 

The rather surprising result that a small variation of 
vertical phase u> can strongly affect the calculation of dy- 
namical friction coefficients K2 and K2 is explained by 
strong cancellation that takes place when one averages 
i • Ai over u>: some terms proportional to StVdist that av- 
erage to zero when Sujdist = can become very large after 
averaging when Sadist 7^ 0. In graphical form the same 
issue has already been illustrated in Figure [T2j Also, 
a careful inspection of our calculation of all other scat- 
tering coefficients including ((Ai) 2 )^^ shows that unlike 

(i • A\) U _ T they are not affected by non-zero Sujdist 7^ 
since they do not suffer from cancellation effects when 
averaged over u>. For these coefficients, our calculation 
neglecting the effects of distant interaction presented in 
Appendix [A"l remains valid. 

Another important point to make here is that although 
there is a large difference between the numerical and an- 
alytical Ki and K2, this discrepancy is not always crit- 
ical for the inclination evolution of planetesimals in the 
regime (|24l) . Indeed, as Figures [4j{5] show, K1/K2 ~ 10 4 
for o e ^5 — 10, which according to equation (|18p implies 
that the gravitational friction term becomes important 
for inclination evolution only if /ii 3> i-e. for the 
velocity evolution of massive bodies driven by their in- 
teraction with low-mass planetesimals. 

The issue of distant interaction and its effect on the 
scattering calculation clearly deserves further work but 
we postpone it for a separate investigation. For now it is 
enough for us to just state that distant interactions serve 



as a plausible explanation for the deviations between the 
numerical and analytical calculations of (i ■ Ai) w , T , K2 
and K%- 

5. VELOCITY EVOLUTION OF PROTOPLANETARY CORES. 

Here we use our results as obtained in previous sec- 
tions to understand the velocity evolution of a sparse 
population of protoplanetary cores in the end of oli- 
garchic phase, when their orbits become crossing — a 
situation described in Jl] For simplicity we will as- 
sume the masses of all cores M c to be roughly equal, i.e. 
Mi = ^2 = A* = Mc/M* and <J{ e ,i},i = °{es},2 = 
which allows us to omit the term proportional to K2 in 
the equation for inclination evolution; see §4 4.41 Also, 
we will neglect the effects of multiple scattering on the 
inclination evolution, which is justified to some extent by 
the small number of bodies involved (and the small num- 
ber of their close encounters); see 8 344.31 We can then use 
our analytical expressions for the scattering coefficients 
to study velocity evolution. 

We assume that initially a e — a e o S> 1 and = o^o 
satisfies constraint ([Ml) . Equation (Tig]) supplemented 
with expressions for the scattering coefficients yields the 
following set of velocity evolution equations: 



dt "° e e M ' 

d(J i ni ,-l ~ 5 

~~dt = ' 
where evolution time T e is given by 



T, 



1 



/' 



nN p R 2 H 



(40) 
(41) 

(42) 



150 yr 



0.01 M ff 



1/3 



AU 

a 



1/2 



30 g cm 



for M* = Mq, and coefficients C e , Ci are constants of 
order unity (their values can be found from equations 
(|33"1) - (13510 . Note that these equations apply equally well 
both to the dense population of planetesimals with over- 
lapping orbits and to the sparse population of cores with 
crossing orbits. 

One can understand the origin of these equations qual- 
itatively based on the fact that scattering coefficients 
are dominated by the large-angle scattering events in 
the thin-disk case, i.e. those that require impact pa- 
rameter Irain ~ Rh/vq. Since we consider the case 
e 3> 1 the relative velocity between planetesimals is 
~ VoVIRh- Then a mean time between encounters of 
a given body with other bodies at an impact parame- 



ter L 



is T, 



zr/2 



Vo/iflNpRjj), where = S p /M c 



is the planetesimal (or core) surface number density. 
When a large-angle scattering event occurs e 2 changes 
by ~ e 2 , while i 2 changes by ~ * 2 ^Q' see ^4 4.31 Then 
dal/dt ~ e 2 /T 7T / 2 , while daf/dt ~ i Vo/T n / 2 , which re- 
sults in equations (f40|) - ((4Tl) if we recall that % ~ a e . 
Integrating equations (f40 |> - (|4T j) we find 



<T e + C e — , 



Oi = a t0 exp 



6C P 



1 I ~6 ~6\ 



(43) 

(44) 
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For KT e one has a e ss a e o and 



<7i = a i0 exp 



s-5 



(45) 



One can see from these solutions that growth of er,; 
has a rather explosive character: <7j increases exponen- 
tially and at the very beginning the inclination growth 
timescale is ~ &eoT e <C T e since a eQ > 1. As a result, 
while <Ji grows by several orders of magnitude, <r e does 
not change that much. This means that protoplanetary 
cores should very rapidly transition from the very thin, 
almost 2D configuration to a vertically extended disk in 
which the condition (f2"4")) is no longer fulfilled. 

In reality inclination will grow not in a continuous fash- 
ion as described by equations (|4"4"|) and (|15"]) but more in a 
step-like way by a factor of ~ Vq as large-angle scattering 
events occur at time intervals of order T v i%. Continuous 
description is going to be useful only after a large num- 
ber of large-angle scattering events have taken place (and 
during this period i would have grown a lot). 

Onset of the core orbit crossings not only increases 
the cores' inclination but also makes possible collisions 
between cores leading to their growth. While the disk 
is geometrically thin, the collision probability of cores 
is rather large, and one may wonder whether the cores 
would grow appreciably during the short period of time 
while the condition is still fulfilled. What matters for 
the core growth during this period is both the accretion 
rate (which is very high) and the time during which the 
inclination of the disk is still in the regime (pM)) , which is 
short. Careful analysis shows that the relative core mass 
increase during the "thin disk" epoch of the core popu- 
lation evolution is very small, much less than unity, pro- 
vided that vo is less than the escape speed from the core 
surface. A simple reason for this is that when the rela- 
tive speed of bodies is less than their surface escape speed 
Imin is larger than the impact parameter leading to colli- 
sions between the cores, so that the inclination strongly 
increases before cores have had a significant chance to 
collide. Thus, all of the late core agglomeration resulting 
in present day terrestrial planets occurs only after the 
disk has become geometrically rather thick, i.e. already 
when Hi ~ a e . 

A similar picture of velocity evolution - relatively 
rapid growth of inclination compared to the growth of 
eccentricity — is expected also in the rather general 
situation of a planetesimal disk transitioning from the 
shear-dominated to the dispersion-dominated dynamical 
regime as a result of gravitational scattering. Indeed, in 
the shear- dominated case eccentricity growth is expected 
to be much faster than the growth of inclination because 
of the geometry of gravitational scattering of planetesi- 
mals in a dynamically cold disk. Thus, when a e becomes 
comparable to unity and continues to grow the condi- 
tion (p2"4"| is fulfilled, meaning that very soon after leaving 
the shear- dominated regime the planetesimal disk should 
rapidly increase its inclination in accordance with equa- 
tions (|4"3"| - (|45p . This qualitative picture has indeed been 
observed in calculations of planetesimal velocity evolu- 
tion based on direct N-body simulations (Ida & Makino 
1992). 

6. DISCUSSION. 



In this paper we have explored for the first time a 
rather special dynamical regime of planetesimal veloc- 
ity evolution represented by the condition (j2"4")l . Previ- 
ous work towards understanding planetesimal dynamics 
has been primarily focused on thick planetesimal disks 
with i ~ e (e.g. Stewart & Ida 2000; Ohtsuki et al. 
2002). Thin disks have been considered by Palmer et 
al. (1993) but their study assumed a razor-thin disk, i.e. 
i = 0, which precluded them from studying a very im- 
portant aspect of the problem - excitation of inclination 
in a thin disk. They have explored horizontal velocity 
excitation in a 2D disk, however, their results cannot 
be directly compared with ours: Palmer et al. computed 
quantities like dv^ jdt - growth rate of the radial velocity 
dispersion of planetesimals — which cannot be directly 
related to our da\jdt since for the latter one also needs 
to know the growth rate of azimuthal velocity disper- 
sion. Nevertheless, their dv^/dt scales linearly with er e , 
in agreement with our equations (l29|) . ([30]) and ([33]) . (f34j) 
for eccentricity-based scattering coefficients. Ida (1990) 
has also recovered a linear dependence of 2D horizontal 
excitation rates on a e numerically. 

The transition between the thin-disk and thick-disk 
regimes of scattering which occurs at i ~ i cr u has not 
been previously investigated. It is known from the stud- 
ies of thick planetesimal disks that scattering coefficients 
tend to diverge as i — > (e.g. Stewart & Ida 2000). On 
the other hand, Ida (1990) has found numerically that in 
a purely 2D disk horizontal scattering coefficients are fi- 
nite, which led him to a conjecture that these coefficients 
should change discontinuously at i = 0. Our present 
study shows this not to be the case. Instead, 3D rates 
increase with decreasing i until i reaches i cr u at which 
point the geometry of scattering changes and scattering 
coefficients smoothly transition to their 2D values. This 
process is best illustrated in Figure [3] where we plot both 
our thin disk results and the scaling of scattering coeffi- 
cients with i in the 3D regime. Using analytical expres- 
sions for various scattering coefficients in Stewart & Ida 
(2000) we have verified that the magnitudes of 3D scat- 
tering coefficients coincide (up to numerical constant of 
order unity) with values of our 2D scattering coefficients 

at l ^ Icrit ■ 

Our final comment concerns multiple scattering en- 
counters resulting in temporary captures such as the or- 
bit shown in Figure [T0l A long time spent by one body 
in the vicinity of the Hill sphere of another opens up the 
possibility of capturing this body into a distant satellite 
orbit if some weak additional perturbation (e.g. gas drag, 
collision with/or gravitational perturbation by an addi- 
tional passing planetesimal) affects the mutual orbit of 
the bodies. By distant satellite we imply a satellite with 
separation larger than the mutual Hill sphere of the two 
bodies, and it seems plausible that the formation of such 
a configuration should somehow involve a high velocity 
encounter between the two objects (relative speed of a 
distant satellite and its parent body exceeds Hill velocity 
CIR H ). 

Multiple scattering encounters between low-velocity 
planetesimals potentially leading to the formation of 
satellites with separations less than Rh have been stud- 
ied by e.g. Iwasaki & Ohtsuki (2007) and Schlichting & 
Sari (2008). The high- velocity regime of multiple scat- 
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tering and temporary capture has not yet been explored 
theoretically and is likely to differ from the low-velocity 
regime. In particular, Schlichting & Sari (2008) found 
that the probability of forming temporary satellite sys- 
tem drops exponentially with the duration of the tem- 
porary capture in the low-velocity regime, and that es- 
sentially no temporary capture systems should exist for 
more than several tens of f2 _1 . However, in the high- 
velocity case we are finding orbits like the one displayed 
in Figure [TU1 which exhibit temporary capture for more 
than 10 3 ft _1 . 

Distant satellites have not yet been discovered in plan- 
etary systems but their dynamics were investigated the- 
oretically by a number of authors (Jackson 1913; Lidov 
& Vashkov'yak 1994a, b). In particular, recently Shen 
& Tremaine (2008) have demonstrated using a map- 
ping approach that distant satellites around some plan- 
ets (Jupiter, Uranus, and Neptune) are stable on time 
scales comparable to the life time of the Solar System. 
Temporary captures resulting from multiple scattering 
of dispersion-dominated planetesimals described in i)44.3l 
present one of the possible ways in which such distant 
satellites may be formed. The dependence of the effi- 
ciency of this formation channel on the dynamical state of 
the planetesimal disk may provide us with an important 
probe of the dynamical characteristics of the early Solar 
System. Needless to say, issues like planetary migration 
or possible chaotic epochs of the dynamical evolution of 
the Solar System planets (Tsiganis et al. 2005; Gomes 
et al. 2005) must significantly complicate the interpre- 
tation of future detection (or non-detection) of distant 
satellites. Nevertheless, the investigation of their forma- 
tion efficiency in temporary capture events like the one 
shown in Figure [10] is a worthwhile exercise. 

7. SUMMARY. 

We investigated the dynamical evolution of vertically 
thin, dispersion-dominated planetesimal disks with ec- 
centricities and inclinations obeying the constraint (|24p . 



In this regime of orbital parameters planetesimals see an 
anisotropic flux of incoming bodies (unlike in the case of 
thick disks), which dramatically changes the character 
of gravitational scattering. In particular, planetesimal 
velocity evolution is dominated by large-angle scattering 
events, unlike in the thick disk case. We derived ana- 
lytical expressions for the scattering coefficients in the 
thin disk regime and compared them with numerical in- 
tegrations of test orbits in the Hill approximation. We 
found good agreement between the two approaches for 
the eccentricity scattering coefficients, while the numer- 
ical inclination scattering coefficients significantly differ 
from their analytical analogs. We demonstrated that this 
discrepancy is caused by the important role of multiple 
scattering events not captured in our analytical calcu- 
lations, and by the distant interactions of planetesimals 
in their approach phase before close encounter. Based 
on these results we have studied the velocity evolution 
of a population of protoplanetary cores in the end of 
the oligarchic phase and shown that the initially small 
inclination of this population grows very rapidly (expo- 
nentially) on a very short timescale. The results of this 
work are useful for understanding the velocity evolution 
of shear-dominated planetesimal disks at the transition 
to the dispersion-dominated regime and for the forma- 
tion of distant satellites of planets. 
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APPENDIX 

SCATTERING COEFFICIENTS IN THE 2D REGIME. 

To compute and analyze scattering coefficients characteristic for thin planetesimal disks we utilize an approach 
developed in Nakazawa et al. (1989), Ida et al. (1993), Tanaka & Ida (1996). For given h, e, and i there are two values 
of the horizontal phase t c ± and time t c ± 



± 



h 2 



1 



1/2 



(h/e) 



±- 



e 2 



3 V h 2 



1 



1/2 



(Al) 



corresponding to the relative orbit passing through the origin (i.e. x = y = 0) in the zero-inclination case (i.e. z 
identically equal to zero). In the case of i ~ e, passage through the origin also implies that z = resulting in a 
constraint on lo. However, in the case of a very thin disk with very small but non-zero inclination all values of u) 
correspond to roughly the same separation from the origin, which is mainly determined by the value of r. Orbits 
passing close to the origin can be expanded about rf- in terms of tj^ — t — with the result that prior to interaction 
one planetesimal approaches another with velocity (scaled in Hill units by QRh) vq = (vo,x, ^o,yy vo yZ ) given by 

1/2 



5* 



0..r 
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5 t,„ = «' h > 



icos(tf — uj), 



and moving on a straight line orbit with the following coordinates of the point of closest approach: 



x± = ± 



-- 1 

vi 
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e 2 



2 k ~ 



i sin(tf — uj), 



where vq = |vq| = [e 2 — (3/4)/! 2 ] 1 ' 2 . Impact parameter of the approach trajectory scaled by Rh is 



~±N21 1/2 
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(A2) 



(A3) 



(A4) 



In all these expressions we have neglected terms higher order in i. 

In the two-body approximation that we adopt here gravitational interaction of planetesimals changes the straight 
line trajectory into a hyperbola defined as 

lcos9 „ 1 

r= . Q , 7> tan0= 7 -2, (A5) 

sm W + cos / i^g 

where 28 is the bending angle of the trajectory (angle between the incoming and outgoing asymptotes of the orbit) 
and / is the true anomaly of the orbit (angle between the line of focii and a particular point on a hyperbola), varying 
from 7r/2 + 9 (incoming) to —n/2 — 8 (outgoing). It is trivial to show that in (x, y, z) = {x\, X2, £3) coordinates this 
hyperbola can be represented as 



x 1 



VQ,; 



■ cos(/ -8)-^- sin(/ - 6>) 



I Vq 

Also, conservation of angular momentum allows one to relate / and t via 

df lv Q 



i = 1,2,3. 



(A6) 



(A7) 



dt r 2 

To compute the changes of orbital elements from equations pU)l - p3|) we also adopt approximation of "instantaneous 
interaction" meaning that we keep time t fixed (and equal to i c ,±) throughout the scattering process. This approx- 
imation works well in high velocity encounters like the ones we are considering here because the interaction time is 
short. It allows us to integrate equation (TT4"|) as follows (and all others in analogous fashion): 

00 ir/2+6 
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2cosC32 n ), Ae y « ^ (oost^? - 2sm^), 



-00 -tt/2-0 

where we cho ose a value of closest 8 to a giv en value of r and then select a value of t c corresponding to t c , see 
equation (|Aip . Then, using equations (|A1|) - (|A8|) one finds that 

(A9) 
(A10) 

(All) 
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8 Ambiguity in the choice of the origin of our r-expansion arises when \r — 
values of r do not produce noticeable contribution to the scattering coefficients. 
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However, trajectories corresponding to these 
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where i = 1,2,3 stands for x, y, z, correspondingly. 

In Figure [8] we compare analytical and numerical results for the changes of various orbital elements in the thin-disk 
limit i « 1 « e, which is of interest for us here. As has been previously shown by Tanaka & Ida (1996), analytical 
results match numerical ones for most values of t if one shifts the origin of r by a small amount rfr <C 1. This shift 
arises from the distant interaction of the two planetesimals before they have experienced close encounter. Such a shift 
in r does not affect in any way our calculation of scattering coefficients averaged over r. 

A striking feature of Figure [8] is the existence of narrow intervals of r in which numerical results strongly deviate in 
seemingly chaotic fashion from the analytical ones. Deviations of Ae, can be of order Ae itself, but the discrepancy 
between the numerical and analytical Ai exceeds analytical value of Ai by orders of magnitude for some values of r. 
The implications of these deviations are discussed in more detail in ^JH 

Neglecting for now this additional complication we average expressions (|A9|) - (|A10[) over u> and r and finally arrive 
at equations (|2l) ]) -([2"g j) , 

ROLE OF DISTANT ENCOUNTERS. 

Let o>o be the initial relative vertical phase of two planetesimals at infinity and to be the value of this phase right 
before the close encounter. Their difference StOdist = lu — u>o is small but nonzero because of the distant interaction 
of planetesimals preceding their close encounter. Previously we assumed to to be equal to ujq (i.e. Su>dist — 0) thus 
neglecting the effec t of di stan t intera ct ion. Let us now see how the fact that StUdist ^ affects calculation of (i • Ai) w . 

Using equations ifAl]) . (fAlj) . ([ATpf . (fATTjl we find 

r A 7~ ~. . ,~ „~o cos(t" — u) cos(i" — Wo) + Vn sin(t" — cj) cos(t" — ujo) /T ^ 1X 

i • Ai = i cos w Az^ + tsmoj Q Ai y = -2i 2 2_ ~ " — (m~2\2 ■ ( Bl ) 

n =± 1 + ^™ w o) 

If we now average this expression over ujo we find (recall that and vo are virtually independent of ujo in the thin-disk 
regime when i -C e) 

(i • Ai) w « -i > ^ - ■ ■„ , (B2) 



Z^ ± 1 + ^)2 



where 
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Sadist = y / -^y^U-i, )<•<*-{/;.' • ^'ohL-u. (B:-!i 



If distant interaction prior to encounter were not taken into account then duidist — and equation (|B2[) would be 
missing the second term in the numerator because averaging over ujq would kill this term completely. However, when 
distant interaction and the possibility of non-zero Sudist are allowed for the omission of the second term may not be 
justified even if j^Wdistl «C 1 because Vq 3> 1 in the situation that we consider. Then it may be possible that the 
product v^SaJdist ^ 1 and dominates the numerator of equation ()B2|) . which makes our neglect of distant interaction 
in calculation of (i • Ai) w unjustified. Thus, distant interaction of planetesimals can indeed explain the discrepancy 
between analytical and numerical values of (i • Ai) w observed in Figure [121 



